Mon. Not. R. Astron. Soc. 000, ITHT31 (20121 Printed 8 October 2012 (MN M£X style file v2.2) 



Simulating Sinking Satellites with Superbox-10 

R. Bien*, T. Brandt, A. Justf 

Astronomisches Rechen- Institut, Zentrum fur Astronomie, University of Heidelberg, Monchhofstr. 12-14, 69120 Heidelberg, Germany 



Accepted 2012 October 3. Received 2012 July 27; in original form 2012 January 23 



ABSTRACT 

Superbox-10 is the successor of Superbox, a particle-mesh code where additional 
grids and sub-grids are applied to regions of high particle density. Previous limitations 
have been solved. For instance, the vertical resolution is improved considerably when 
flattened grids are used. Since the computationally most intensive part is the Fast 
Fourier Transform, we introduce a parallelised version using the library fftw, result- 
ing in a speed-up of a few. The new features are tested using a galaxy model consisting 
of an exponential disc, a bulge and a dark matter halo. We demonstrate that the use 
of flattened grids efficiently reduces numerical heating. We simulate the merging of 
disc-bulge-halo galaxies with small spherical satellites. As a result, satellites on orbits 
with both low eccentricity and inclination heat the disc most efficiently. Moreover, we 
find that most of the satellite's energy and angular momentum is transfered to the 
halo. 

Key words: methods: numerical - stellar dynamics - galaxies: evolution - galaxies: 
kinematics and dynamics - galaxies: structure. 



1 INTRODUCTION 

Observations of stars in the solar neighbourhood show 
an increase of the velocity dispersion (the random mo- 
tion) proportional to f - 3 - - - 6 ^ where t is the stellar age 
(Wielen 1977; Holmberg, Nordstrom & Andersen 2007). 
This effect can be interpreted as a heating of the Galactic 
disc with time. However, it is unclear what is the cause 
of this heating. Several mechanisms have been proposed: 
Massive black holes in the gala ctic halo as a source o f 
heating were investi gated by lLacev fe Ostriker 
Wielen fc Fuchs (19901), iRix fc Lake ( 19991 ) and 



Hanninen fe Flvnn (2003 ) , but seem to be excluded 



by observations. lArdi. Tsuchiva fc Burkert (20031 ) 
considered massive clumps of dark matter as a 
possible heating agent. Giant molecular clouds 
were found to cause mos tly vertical heating b y 
Lacev (198 1) and late r iHimninen fc Flvnn (2003] . 



Carlberg fc Sellwood (19851 ) and ICarlberg (19871 ) show 



that transient s piral waves heat efficiently in the G alactic 
plane. Lastly, Quinn. Hernquist fc Fullagar (19*931 ) and 



iVelazquez fc White (199g T simulate merging with small 
satellites and find both radial and vertical heating. This 
last mechanism is discussed in this paper. It should be 
mentioned, however, that radial mixing may influence the 
age-velocity dispersion relation significantly (Schonrich & 
Binney 2009; Sellwood & Binney 2001; Roskar et al. 2008). 



* E-mail: reinhold@ari.uni-hcidelberg.de 
f E-mail: just@ari.uni-hcidclberg.de 



The merging of galaxies can be divided into three 
categories according to their mass ratios. Major merg- 
ers with mass ratios 1:1 to 1:3 of the total galaxy 
masses usually destroy discs and result in an early 
type remnant (|Naab. Jesseit fc Burkert 20061 . and ref- 
erences therein). Minor mergers with mass ratios 1:3 
to 1:1 of the total galaxy masses d estroy a thin 
disc (|Purcell. Bullock fc Kazantzidis 201~bT ). The minor 
merger events are discussed in the context of the for 



mation of thick discs JVillalobos fc Helmi 200sll 200 
and substructure in the halo ( Newberg et al. 200 
IPurcell. Bullock fc Kazantzidis 2010l ). The third category 
is the merging of satellite galaxies (or dark matter clumps) 
with smaller mass ratios. These merger or interaction events 
can be characterized as a perturbation of the primary galaxy 
and it is more convenient to quantify the mass ratio in 
units of the disc mass of the primary. In this context the 
survival of a thin disc over 10 Gyr in a ACDM universe lead 
to serious constraints on the d umpi ness of the DM halo 
jArdi. Tsuchiva fc Burkert 20031 ). iHopkins et al. (2008! ) 
compared the efficiency of the disc by satellite galaxy 
mergers of different authors and argued that it scales with 
the square of the satellite mass. This is consistent with the 
interpretation heating by dynamical friction during each 
disc crossing event. As a consequence, the upper mass end of 
the perturbers are most efficient and determine the survival 
of a thin disc. Kazantzidis et al. (2008, 2009) performed 
a comprehensive study of the impact of DM subhaloes on 
a Milky Way like disc galaxy in a realistic cosmological 
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context. They investigated the resulting substructures and 
kinematic features in disc and halo. 

In numerical simulations the long-term (or secular) 
dynamical evolution of a thin galactic disc is very sen- 
sitive to the spatial resolution, the level of numerical 
noise, and the dynamical feedback with the DM halo. 
A live dark matter halo is important for two reasons. 
IVelazouez fc White (HH ) have shown that the reaction 
of the disc on the perturbation of a satellite galaxy is 
significantly higher with a live halo compared to a rigid 
dark matter potential. Secondly, dynamical friction in the 
dark matter is important for the orbital evolution and 
thus energy loss of satellites with masses exceeding ~ 
10 8 Mq. This is exactly the mass range dominating the 
heating rate of the disc by satellite galaxies. Addition- 
ally, the mass of the dark matter particle must not be 
too large in order to avoid numerical heating of the disc. 
In most simulations the disc is represented by less than 
1 million particles leading to a significant thickening of 
the disc by numerical two-body relaxation, which is much 
stronger in direct A-bod y codes and Tree codes co mpared 
to particle mesh-codes. I n lVelazquez fc White (19991 ) and in 
lHavashi fc Chiba (20061 ) the nu merical heating of the unper - 
turbed disc was shown explicitly. I Velazquez fc White (19991 ) 
used a differential method to quantify the dynamical heat- 
ing of infalling satellites. This may underestimate the heat- 
ing rate, since a more realistic thinner, dynamically cooler, 
disc is more sensitive to perturbations. On the other hand, 
the feedback on the satellite galaxy may be larger by the 
stronger disc shocking event. Particle-mesh codes are pre- 
destined for the simulation of collisionless systems such as 
galaxies, because two-body relaxation is strongly suppressed 
by the partial decoupling of the point-like structure in or- 
bital motion of t he particles and the gravi tational potential 
based on the grid. lKhoperskov et al. (20071 ) have shown that 
the growth of global mode perturbations are easily damped 
out by noise due to an insufficient particle number. If the 
excitation of spiral structure or bar-like perturbations are 
an important mechanism for the energy transfer from the 
satellite to the disc, then a very high resolution well above 
1 million particles is needed to reproduce a realistic heating 
rate. 

Improving the resolution of discs in vertical direction 
has always been a challenge. There is, indeed, a further prob- 
lem which is inherent in any numerical technique that is used 
for the simulation of disc galaxies. Stellar discs have relax- 
ation times larger than a Hubble time, i. e. they are collision- 
less systems. This implies that in simulations the thickness 
of an unperturbed disc indicates how well the code would 
model a collisionless gravitating system. The term "unper- 
turbed" denotes the absence of perturbations such as stellar 
bars, spiral arms, or molecular clouds. One often notes, in 
simulations, that discs become thicker with time. This is 
caused by the graininess of the particle distribution and the 
limited spatial resolution. We refer to this effect as "numer- 
ical heating". 

For two decades Superbox is used as a successful tool 
to simulate the dynamics of isolated and interacting galax- 
ies. A first description was given by Bien, Fuchs & Wielen 
(1991), together with a consideration of the direct A-body 
technique and a tree- code. For a more detaile d discussion the 
reader is referred to IFellhauer et al. (20001 ). We also men- 




Figure 1. Projection of a simplified galaxy model in a 3D box. 
Shown are the basic parameters -Rcore, Rout and -R sys . 



tion Fellhauer's lecture notes (|Fellhauer 20081 ). The code was 
evol ving from the conventional p article-mesh technique (see, 
e. g.. lHocknev fc Eastwood 19881 ) by applying grids and sub- 
grids in regions of high particle density. The advantage of 
Superbox is that the code can run on any workstation or 
PC, giving reliable results. The code thus proves to be a se- 
rious alternat ive to direct A-body sch emes and tree-codes. 

Initially, iMadeiskv fc Bien (1993h applied Superbox 
to compare numerical experiments to observations of the 
high-velocity encounter of NGC 4782/4783. They found 
a convincing description of the morphological and kine- 
matical structure, and estimated the time elapsed since 
the closest approach of the two interacting galaxies. This 
was the beginning of a series of research projects using 
Superbox as a tool. We note as examples the dynamica l 
evolution of a satellite galaxy (|Klessen fc Kroupa 19981 ) . 
the decay of dwarf galaxies in dark matter haloes 
(|Penarrubia. Kroupa fc Boilv 2002 | ). and dynamica l fric- 



tion, compare Penarrubia. Just fc Kroupa (20041 ) and 



Just fc Penarrubia (20051 ). It should be mentioned that 
Spinnato et al. (20031 ) studied t he inspiral of a black 



hole to the Galactic centre, and iKhoperskov et al. (2007 ) 
considered unstable modes in a disc. IFellhauer et al. (2008 ) 
published on the dynamics of the Bootes dwarf galaxy. 
IPefiarrubia. Navarro fc McConnachie (20 08) investigated 
the tidal evolution of dwarfs of the Local Group. Most 
recently, a study was made on the dy namical friction o f 
massive objects in galactic centres by ljust et al. (20lJ ). 
Our list is by far not complete. 

So far, Superbox was not able to adapt the sub- 
grids to discs. A higher resolution in vertical direc- 
tion can be ach ieved by introducing flattened grids, see 
iBien et al. (20081 ). This new feature and the application to 
dynamical heating of galactic discs are the main topics of 
the present paper. The improved code is called Superbox- 
10. As Superbox, Superbox-10 is rather free of numerical 
heating. 

Below we first describe Superbox in detail and then 
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extend to Superbox-10. It follows a discussion of the par- 
allelised version which results in a significant speed-up. Then 
our test of improved vertical resolution is described. As an 
important application we consider disc heating by satellite 
galaxies. In particular, we address the ratio of radial and 
vertical velocity dispersion and the transfer of energy and 
angular momentum. Special attention is given to the inter- 
action of the infalling satellite with the dark matter halo. 

In a conference proceedings publication 
l|Bien et al. 20081 ) we presented preliminary results on 
an improved version of Superbox with increased vertical 
resolution of the disc: Negligible disc thickening of an 
isolated disc after 2.5 Gyr. The present paper should be 
seen as a pilot study for high resolution simulations of disc 
heating by satellite galaxy mergers. Two new aspects are 
presented: (1) No correction for numerical heating is needed 
and we can measure the heating rate directly. Dynamical 
friction in the dark matter halo is included automatically. 
(2) We analyse for the first time additionally the transfer 
of angular momentum. This is important for understanding 
the physical processes responsible for the disc heating. For 
this pilot study we select a satellite with a relatively high 
mass on eccentric orbits, which is destroyed after a few disc 
crossing events. 
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Figure 2. The five grids of Superbox. Grid 1 has the highest 
resolution and resolves the core. Grid 1, grid 2 and grid 3 move 
through the local universe, defined by grid 4 and grid 5. 



2 A REVIEW OF SUPERBOX 
2.1 Basic concepts 

The conventional particle-mesh technique considers a set 
of massive particles, often called "superstars", in a three- 
dimensional Cartesian grid ( "box" ) consisting of TV x TV x TV 
cells which represents a relevant part of the universe. In this 
section, the length of each cell is the unity. We suppose that 
p a ,b,c is the mean density in the cell with indices a, b, c which 
are integers running from to TV — 1. Then the potential can 
be obtained by solving Poisson's equation 

JV-l JV-l JV-l 

$a,b,c = G ^ 2J ^ P(,t,X H i-a,r,-b,(-c (1) 
{=0 77=0 C=0 

Here, H(- atV -bx-c 0- e - Green's function) is the inverse dis- 
tance between the points having indices a, b, c and f, 77, £, 
namely 

ff e - a ,„_ b ,C-c = [« - a) 2 + fa - bf + (C - c) 2 ]- 1/2 (2) 

We note that the term Gp^ iV ^H^- atV -i,^- c accounts for the 
contribution of the cell f , 77, £ to the cell a, b, c. Thus 4v,t.,c 
gives the correct potential assigned to cell a, b, c. The evalu- 
ation of all TV x TV x TV values <& a ,b,c takes a time proportional 
to the number of cells squared, (TV x TV x TV) 2 . 

The direct convolution can be replaced by discrete 
Fourier transforms where TV is supposed to be a power of 
2. First p a ,b,c and H a ,b,c are transformed, resulting in 

N-1N-1N-1 ^ 

Pk,l,m = ^ ^ ^2 Pa,b,cexp[-i — (ka + lb + mc)] (3) 

a = 6=0 c=0 

and Hk,i,m, analogously. The potential is then obtained by 
the inverse transformation of the product Pk,i,m,Hk,i,mt 

$ ,6,c = 



N-lN-lN-l ^ 

jp ^ ^ ^ Pk,l,mHk,i, m exp[+i — (ka + lb + mc)} (4) 

fe = i = m=0 

The procedure can dramatically be accelerated when a Fast 
Fourier Transform (FFT) is applied. The computing time is 
then roughly proportional to TV 3 log TV. 

The potential, as derived so far by Fourier transforms, 
is correct for periodic systems only. The exact potential of 
isolated systems can be obtained by doubling the number of 
cells in each coordinate direction, i. e. a grid is considered 
which contains 2TV x 2T V x 2TV cells. A rigorous proof can 
be found in the paper bv lEastwood fc Brownrigg (19791 ). As 
before, the TV x TV x TV sub-grid contains the particles ( "active 
region"). The remaining space is left empty. Likewise, H is 
extended over the 2TV x 2TV x 2TV cells such that it satisfies 

H2N-a,b,c = H2N-a,2N-b,c 
= H2N-a,b,2N-c = H2N -a ,2N -b,2N - c = H a ,2N-b,c 

= H a ^2N-b,2N-c = Ha,b,2N-c = -ffa,b,c (5) 

for ^ a, b, c ^ TV. Having transformed the extended func- 
tions p and H, the correct potential $ can be found in the 
active grid. Outside, $ is unphysical. 

Numerical differentiation of the potential with respect 
to the coordinates leads for each particle j to the accelera- 
tion (a x )j, (a y )j, (a-z)j- The leap-frog scheme is applied to 
integrate the equations of motions. In ^-direction the algo- 
rithm reads 

xj(t + At/2) = ±j(t - At/2) + At(a x )j 

Xj(t + At) = Xj(t) + AtXj(t + At/2) (6) 

where t is the time and At denotes the constant integration 
step. 

The particle-mesh technique holds some subtle twists. 
For instance for H a ,b,c it suffices a(TV+l)x(TV+l)x(TV+l) 
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grid. From the symmetry conditions ([5]) follows that the 
function is then known on the whole 2N x 2N x 2N grid. The 
sub-grid is overwritten by H which shows similar symme- 
tries. H is calculated only once, since it is fixed over the in- 
tegration period. For each integration step the grid contain- 
ing the density p is overwritten by the Fourier transform (p 
which, in turn, is overwritten by pH). This product is finally 
overwritten by The advantages of the method are obvi- 
ous. The conventional particle-mesh technique was orders of 
magnitude faster than direct iV-body codes at the time. This 
allows the integration of a large number of particles. Thus, 
statistical noise is extremely low. Particles leaving the out- 
ermost grid are lost. The user should thus carefully consider 
the number of particles involved in the computation. 

The greatest disadvantage, however, is the low spatial 
resolution in regions of higher particle densities, e. g. the 
cores of spherical galaxies. An improvement is Superbox 
which treats the particles self-consistently in a nested sys- 
tem of grids. Here, the term "self-consistent" means that all 
particles have masses which define the total potential. For 
simplicity, we consider a spherical galaxy of radius _R ou t ■ Fig. 
1 shows how the particles are subdivided into a central re- 
gion ( "core" ) where the distance r from the centre is < R COI<! 
and a region with R C orc < r < R ou t- The local universe. is 
given by it sys . 

As shown in Fig. 2 Superbox utilises five grids. They 
all have the same number of cells, TV 3 . 

(i) Grid 1 has the highest resolution and resolves the core 
which contains all particles with < -Rcore- 

(ii) Grid 2 contains the core, too, but the resolution is 
intermediate. 

(iii) Grid 3 is equal to grid 2, but contains only particles 

With J?core < r < -Rout ■ 

(iv) Grid 4 is the fixed global grid. It contains only par- 
ticles inside R ou t- 

(v) Grid 5 is equal to grid 4 with all particles outside 

Rout ■ 



where [ ] denotes the nearest integer function. The factor 7 
enhances (or shrinks) an area. If 7 = 10 is assigned to grid 
1, the resolution increases by a factor of 10. Inside the cell 
with indices i, j, k, let the particle's coordinates be Ax, Ay, 
and Az. Then its acceleration, e. g. in ^-direction is found 
by the expression 



t - $i_i 



+ 



2l x 

$»i+i,j,fc + $i-i,j,fc - 2$i,j, fc . 

72 ^ X 
1 X 

&i+l,j + l,k — $i—l,j+l,k + &i—l,j—l,h — ^i+l,j-l,fc a 

47 7 y 

"-tlx <■■ y 

$i+i,j,fc+i — 3>i-i,j,fc+i + $i-i,j,fc-i — $j + i,j,fc_i ^ , . 



where l x — l y = l z is the length of the cell. This scheme re- 
quires two empty cells at each boundary, i. e., only N— 4 cells 
per dimension are considered. We note that in Superbox, 
the three lengths are equal. When all velocities are updated 
by applying these accelerations, then new positions of the 
particles are calculated. After that, a new integration cycle 
starts. 



2.2 Extensions to Superbox-10 

For version 10, Superbox's original Fortran 77 code has 
been ported to Fortran 95 and arranged into several mod- 
ules, to make future extension easier. 



2.2.1 Individual masses 

The conventional particle-mesh technique considers only a 
single mass for all particles. In Superbox- 10, an individual 
mass can be assigned to each particle. In practice, however, 
the particles of each subgroup (i.e. halo, disc, etc.) carry 
identical mass. 



As the galaxy moves, grid 1, grid 2 and grid 3 move through 
the local universe (i. e. grid 4 and grid 5, respectively). This 
is done by centreing the inner and intermediate grids on the 
density maximum of the galaxy in question. Alternatively, 
the centre of mass can be considered. Particles outside R sys 
are lost, compare the discussion above. The code is not re- 
stricted to one galaxy. Since the potentials, and thus the 
accelerations, are additive, all galaxies are treated sequen- 
tially in the same five grids. The corresponding five total 
potentials are $1, $2, $3, $4, and $5 

(i) For a particle with r < R CO re, the correct potential is 
*i + $3 + $5- 

(ii) For a particle with, 7? C oro < r < Rout the sum $2 + 
$3 + $5 is taken. 

(iii) When the particle is outside i? ou t then $4 and $5 
are used. 

In principle, number and type of the galaxies is arbitrary. 

Two features should be emphasized. A particle at the 
point (a;, y, z) is assigned to the cell with indices 



(7) 



2.2.2 Improved vertical resolution 

Superbox poorly resolves stellar discs in vertical direction. 
Superbox- 10 overcomes this shortcoming. We explain the 
basic idea by taking the example of a disc-bulge-halo galaxy. 
First, the potential of the disc-bulge-halo galaxy is calcu- 
lated. The intermediate grids grid 2 and grid 3 are flattened 
along the corresponding z-axis. When, for instance, the flat- 
tening is q = 1/4 the resolution is improved by a factor of 4. 
As a consequence, equations [5] [7] and [5] need to be changed. 
Yet, there is a restriction to q. In order to cover grid 1 and 
at least two cells of grid 4 (and 5, respectively), q should not 
be smaller than 

Rcorc 4 Rs 



9crit 



(Rcc 

VRo 



ut ' N — 4 Rout 



(9) 
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r AH 
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lx + -_ 
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_7V + y_ 


, c — 





Otherwise, spurious results are obtained. If more than one 
gala xy is considered, q < 1 applies only to the first galaxy. 

iBien et al. f200Sfl simulated a disc-bulge halo galaxy 
using Superbox-10. The authors made a compaxi- 
son with a code based on the TREE-GRAPE scheme 
|Fukushige. Makino fc Kawai 20051 ). For 2,577,235 particles 
the CPU time turned out to be about three days on a single 
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Table 1. Maximum and achieved speed-up, S ma x and S, as func- 
tion of grid size N and particle number n. f = tFFT/itoti see 
text. 



TV 


n 


/ 




s 


S/ Smax. 


64 


500k 


0.62 


2.61 


2.48 


0.95 


128 


500k 


0.92 


12.58 


8.96 


0.71 




2M 


0.74 


3.80 


3.26 


0.86 


256 


500k 


0.99 


71.56 


23.49 


0.33 




2M 


0.95 


20.20 


11.90 


0.59 




10M 


0.78 


4.48 


4.05 


0.90 



GRAPE6a board. The CPU time of Superbox-10 (which 
is still called Superbox in that paper) on a customary PC 
is of the same ord er. This result sho ws impressively the effi- 
ciency of the code. iBien et al. (2008) came to the conclusion 
that the extended code 

(i) is very fast, 

(ii) does not produce noticeable numerical disc heating, 
and 

(iii) allows an improved vertical resolution. 
2.2.3 Parallelisation 

The computationally most intensive part of Superbox-10 is 
the calculation of the potential using the Fast Fourier Trans- 
form. The non-parallelised version of Superbox- 10 applies 
the 1D-FFT routi ne realft, taken from Numerical Recipes 
(|Press et al. 19921 ). The 3D-FFT is calculated by doing TV x 
TV ID transforms in each direction. In the parallelised ver- 
sion, this scheme has been replace d by a 3D-FFT from th e 
fftw library (version 2.1.5), see iFrigo fc Johnson f2005h . 
This routine divides the 3D array into slices along one di- 
mension and distributes these slices among multiple pro- 
cessors. The processors then jointly calculate the 3D-FFT, 
using the Message Passing Interface to communicate. The 
modular design of Superbox- 10 allows us to replace the 
FFT routine without changing much of the rest of the code. 
The non-parallelised FFT is still available. 

In the following we introduce the fraction / = £fft /ttot 
in order to analyse the non-parallelised version in more de- 
tail. Here £fft is the time spend in the FFT routine per 
integration step and ttot is the total time used for the in - 
tegration step. According to Amdahl's law (I Amdahl 1967t ). 
the maximum achievable speed-up by improving the FFT is 
then given by 

Smax = Y~f ( 10 ) 

Table Q] lists /, Smax and the speed-up S achieved in bench- 
marks. 

The values depend on the grid size TV as well as on the 
particle number n. For n > TV 3 , the integration of the parti- 
cles' orbits starts to dominate over the potential calculation. 
Thus, the higher TV, the greater the expected speed-up will 
be. 

Benchmarks were run on the Titan cluster at the As- 
tronomisches Rechen-Institut, which has 32 nodes, each con- 
taining a Xeon 3.2 Ghz dual core CPU and a nVidia GeForce 
9800 GTX Graphics Processing Unit (GPU). However, we 




1 4 8 12 16 20 24 28 



number of processes 

Figure 3. Speed-up S compared to non-parallelised Superbox as 
function of the number of processors for varying TV 3 , the number 
of grid cells, and particle number n. 

used only one core per CPU in order to have more memory 
available, allowing for higher grid sizes. The GPUs are cur- 
rently not used by Superbox-10, but a version suitable for 
GPUs is in preparation. 

Figure f3] shows the measured speed-up S as function 
of the number of processes in comparison to the non- 
parallelised version of Superbox with the old FFT when 
varying both the grid cell number TV and the particle num- 
ber n. The application of FFTW alone already results in 
a speed-up between 1.8 and 2.4. In cases where the curves 
reach their maximum (TV = 256 for n = 10 x 10 6 and TV = 
128 for all n) the speed-up lies between 0.7 and 0.9 S ma x- 
In the remaining two cases (TV = 256 for n = 0.5 x 10 6 and 
n = 2 x 10 6 ) it lies between 0.3 and 0.6 Smax- When more 
processors are applied one can expect a further increase. 

The fact that the speed-up almost reaches Smax shows 
that there is only a small overhead due to communication 
between the nodes. Thus, from the values of Smax in Ta- 
ble [T] we can formulate the following rule of thumb: When 
using p processors and a grid of size TV 3 , the number of par- 
ticles should not be greater than 4TV 3 /p. Conversely, when 
simulating n particles with a grid of size TV 3 , no more than 
4TV 3 /n processors should be used. 



3 MODELS 

In this section we will describe the models used in our sim- 
ulations of isolated disc-bulge-halo galaxies (section 3} and 
merging of galaxies with satellites (section [5]l . 

3.1 Galaxy model 

Our galaxy model consists of a disc, a bulge and a dark 
matter halo component. The density profile of the disc is 
exponential in radial and isothermal in vertical direction for 

R < -Rmax and \z\ < z maj : 



6 R. Bien, T. Brandt, A. Just 



Table 2. The three components of our galaxy model and the 
satellite. Shown are the number of particles n, the total mass M 
and the mass per particle m. 





n [10 6 ] 


M [10 10 M ] 


m [10 3 M ] 


disc 


5.19 


4.82 


9.30 


bulge 


2.30 


2.14 


9.30 


halo 


2.15 


20.0 


93.0 


satellite 


0.50 


0.54 


10.08 



g(R,z) — go e R ^ h sech 2 (z/.zo) 



(11) 



The disc has a scale length h = 2.5 kpc and a thick- 
ness zo = 0.6 kpc. The profile is cut off at i? ma x = 10 ft 
and 2 max = 10 Zo- This excludes 0.05 per cent of the mass 
an infinite profile would have. At R = 8 kpc the Toomre 
parameter has a value of Q = 2. 

Both bulge and halo have a cropped Hernquist profile 
|Hernquist 1990h 



e(r) = Qo 



r(r + a) 3 



(12) 



with scale radius a and cutoff radius r c . 

The scale radius of the bulge is a = 0.5 kpc and is cut 
off at r c — 14 a. The halo's scale radius is a — 16.8 kpc and 
it is cut off at r c — 5 a. 

We emphasize that both bulge and halo are represented 
by live particles. The importance of treating the dark mat- 
ter halo of galaxies as a live component, i.e., dynamically 
evolving and in mutual interaction with the baryonic com- 
ponent, has been stressed in self-consistent iV-body simu- 
lations of isolated disc galaxies (e.g., Athanassoula 2002; 
Dubinski, Berentzen & Shlosman 2009). The angular mo- 
mentum exchange between the disc and dark matter halo in 
such cases is mediated by dynamical resonances. These pro- 
cesses cannot be resolved when using a static, e.g., analytic 
prescription for the halo potential. 

The m odel has been implemented with the progra mme 
MaGaLie (|Boilv. Kroupa fc Penarrubia-Garrido 20011 ). We 
extended MaGaLie to allow for up to 10 million particles. 

The particle numbers and masses of the three compo- 
nents are listed in Tabled In order to reduce the total num- 
ber of particles, the particles of the halo are chosen to be ten 
times as massive as those of disc and bulge This does not 
pose a problem though, because they are still light enough 
not to cause numerical heating. 

Initially, the model is not completely in equilibrium. 
Therefore, it is evolved in isolation with high resolution 
(N — 256, q = 0.25) until it reaches equilibrium. This resolu- 
tion is high enough not to introduce any numerical heating. 

Figure [4] shows the disc's mean rotational velocity 
and its velocity dispersions or and a z in radial and vertical 
direction after this initial relaxation. Also shown is the cir- 
cular velocity v c . The stellar disc has a flat rotation curve 
and an exponentially decreasing velocity dispersion. 



3.2 Satellite model 

The satellite for the merger simulations has a Plummer pro- 
file 



250 



200 



150 



100 




Figure 4. The disc's mean rotational velocity (v&), radial ve- 
locity dispersion ctr, and vertical velocity dispersion <r z . Also 
shown is the circular velocity v c . Note that the galaxy has al- 
ready reached equilibrium. 



Table 3. Vertical resolution in parsec for various combinations 
of q and N. Note that Rout = 28 kpc. The radial resolution is 
independent of q and has the same value as the vertical resolution 
for 5 = 1. 



q\N 


64 


128 


256 


512 


1 


933 


452 


222 


110 


1/2 


467 


226 


111 


55 


1/4 


233 


113 


56 


28 


1/8 




56 


27 


11 



p(r) = 



3M 
47ra 3 



1 + 



(13) 



with mass M = 5.4 x 10 Mq = 0.11 x M d i ac and scale 
radius a = 1.5 kpc. Its profile is cut off at r c = 10 a. It 
consists of 5 x 10 5 particles (see Table [2} • Like the galaxy, 
it is first evolved in isolation. 



4 ISOLATED GALAXY MODELS 



iBien et al. (2008h showed that Superbox intrinsically has a 
very low level of numerical heating. We discuss now in de- 
tail the dependence of the numerical heating on the vertical 
flattening factor q. To that end we simulate the previously 
discussed disc-bulge-halo galaxy model in isolation with var- 
ious resolutions and measure the change in vertical thickness 
of the disc component. 

4.1 Simulation parameters 

In all simulations, the grids have i? C ore = 3.5 kpc, R ou t = 
28 kpc and R sys = 105 kpc. Depending on the number of 
cells iV 3 and flattening parameter q, the simulations have 
vertical resolutions listed in Table [3] The radial resolution 
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Figure 5. Standard deviation of the z-coordinate of disc particles 
as function of R for varying N and q at t = 1 Gyr, together with 
the initial curve. The values are calculated in radial bins of equal 
particle numbers. 
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Figure 6. Standard deviation of the z-coordinate of disc particles 
as function of R for varying N and q at t = 1 Gyr, together with 
the initial curve. The values are calculated in radial bins of equal 
particle numbers. 



for a certain N = No has the same value as the vertical 
resolution for N = No and q = 1. 

We run low-resolution simulations with N = 64 and 
q — 1, 0.5, 0.25. The radial resolution is 933 pc and the 
vertical resolution amounts to 933, 467 and 233 pc, re- 
spectively. These simulations are compared to a medium- 
resolution simulation (N = 128, q = 0.25) where the radial 
resolution is 452 pc and the vertical resolution is 113 pc. 
Additionally, a comparison to the initial, relaxed, system is 
made. A time step of 0.4 Myr is used and the length of the 
integration corresponds to 1 Gyr. 



4.2 Results 

Figure [5] shows the root mean square (RMS) of the z- 
coordinate of all disc particles, yj (z 2 ), as a function of radial 
distance R. We used radial bins of equal particle numbers. 
The RMS can be taken as a measure of the thickness of the 
disc. The low-resolution simulation (N = 64, q = 1) over- 
estimates the thickness significantly. The maximum devia- 
tion from the initial thickness is about 16 per cent. Flat- 
tening the intermediate grid by a factor of 4, brings the 
thickness down to within about 4 per cent of the initial val- 
ues, without increasing the computation time. The medium- 
resolution simulation (N — 128, q = 0.25) deviates at most 
m 1.6 per cent from the initial values, but is closer than 1 
per cent for most of the radial range. 

To demonstrate the effect of flattening, Figure[6]directly 
compares two simulations with the same vertical resolution 
- N = 128 with q = 0.25 and N = 256 with q = 0.25 - 
corresponding to a vertical cell length of about 112 pc. As 
a reference the case of N = 128 without flattening is also 
shown. As can be seen, introducing flattening in the case of 
medium resolution diminishes thickness, down to about the 
same value as in the case of high resolution. 

Figure [7] displays the vertical density profile averaged 



over the radial interval 7.5 kpc ^ R ^ 8.5 kpc. The initial 
profile (solid line) and the profile of the medium-resolution 
simulation (dash-dotted line) are almost identical. In the 
case of the low-resolution simulations (dashed line) , the pro- 
file is significantly widened for \z\ > 1 kpc (see Figure [7] (b)). 
This is a sign of numerical heating. Reducing q brings the 
profile closer to the initial one. For q = 0.25, deviations are 
only visible beyond 2 kpc. 

Figure [8] shows the time evolution of the RMS for vary- 
ing grid size N and flattening factor q in the region around 
R = 8 kpc. In the first 200 Myr, the model adapts to the 
new grid structure, which causes a change in the RMS. Af- 
ter that, the RMS increases slightly with time in the case 
of the low-resolution simulations. It reaches values that are 
8.7 per cent (q = 1), 3.6 per cent (q — 0.5) and 2.4 per cent 
(q = 0.25) greater than the initial value. In the medium- 
resolution simulation however, it remains stable at a devia- 
tion of about 0.9 per cent. 



5 MERGER SIMULATIONS 

One application that benefits from the increased z-resolution 
is the study of the dynamical heating of galactic discs caused 
by the merging with satellite galaxies. We simulate the merg- 
ing of a small satellite with a disc-bulge-halo galaxy for vary- 
ing initial positions and velocities of the satellite. 

5.1 Simulation parameters 

All satellites are initially at a distance of _Ra = 25 kpc 
from the centre of the galaxy and have a velocity of either 
115 km s _1 or 81.3 km s _1 . If all the mass of the galaxy in- 
side Ra were concentrated in a point mass, these velocities 
would result in elliptical orbits with apocentre distance _Ra 
and eccentricities e = 0.56 and e = 0.78, respectively. Here, 
the eccentricity is defined as 
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Figure 7. Volume density p as function of z (averaged over the 
radial range 7.5 kpc ^ R ^ 8.5 kpc) for varying N and q at 
t = 1 Gyr. The solid line shows the initial values, (a) Whole 
range, (b) Right tail. 



Table 4. Satellite parameters e (eccentricity), v (velocity), Ra 
(apocentre), Rp (pericentre) and i (orbital inclination). 



e 


v [km/s] 


i?A [kpc] 


R P [kpc] 


i [°] 


0.56 


115 


25 


7.05 


0, 10, . . ., 180 


0.78 


81.3 


25 


3.08 


0, 10, . . ., 120 



6 2 



with semi-major axis a and semi-minor axis b. In the two- 
body problem, the pericentre distance of these orbits would 
be Rp m 7 kpc and Rp w 3 kpc, respectively. 

In addition to the eccentricity, we vary the orbital in- 
clination i relative to the plane of the disc between 0° and 
180° for e = 0.56 in steps of 10°, and likewise between 0° 
and 120° for e = 0.78. For i < 90° the satellite is on a pro- 
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Figure 8. Standard deviation of the z-coordinate of disc particles 
as function of time for varying grid size N and flattening q. The 
standard deviation is calculated in a region around R = 8 kpc. 
The solid line shows the initial value. 



grade orbit as compared to the galactic rotation, while for 
i > 90° the orbit is retrograde. For i = 90° the satellite's 
initial velocity is perpendicular to the bulk motion of the 
disc. Table |3] summarises these parameters. 

As shown in the previous section, a grid cell number of 
iV = 128 with flattening factor q = 0.25 causes no significant 
numerical heating. We adopt these values for our merger 
simulations. All simulations run for 1 Gyr. After that, the 
satellites are completely dissolved. As a control, the galaxy 
model is also evolved in isolation. 

At the end of the simulation, the plane of the disc is 
tilted by a few degrees. In the following, the data are evalu- 
ated in a coordinate system where x and y define the galactic 
plane and z is perpendicular to it. 



5.2 Vertical profile 

Figure [5] shows the vertical density profile of the disc after 
merging. Only satellites with significant effect on the disc 
are shown (e = 0.56, i = 0° . . .30°). Compared to the iso- 
lated simulation (i.e. without satellite), the central density is 
decreased by a factor of 2 at the most (i — 20°). In the inner 
part (|z| < 1 kpc), the original sech 2 (z/zo) profile remains 
a good fit, albeit with a greater thickness zq. The isolated 
galaxy has zq = 0.6 kpc, while for i = 20° the thickness is 
increased to 0.8 kpc. In the outer part, however, the density 
is increased as compared to a sech 2 (z/,zo) law. While the 
thickness of 0.8 kpc corresponds to a scale height of 0.4 kpc, 
the outer part is better fitted by an exponential function 
wit h scale height ~ 0.65 kpc . 

iHavashi fc Chiba (20061 ') provide a formula for the in- 
crease of the thickness of the disc caused by cold dark mat- 
ter subhaloes, which should also be applicable to satellite 
mergers. In our notation their formula is given by 



Az ~ 8h 



/ Msat V 

VMd.sJ 



0.25 kpc 
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Figure 9. Volume density of the disc in the region 7.5 kpc 
R ^ 8.5 kpc after the merging (t = 1 Gyr) as function of z for 
different satellite orbital inclinations i. All satellites show have 
initial orbital eccentricity e = 0.56. The solid line corresponds to 
the isolated simulation, i.e., without satellite, (a) Whole range, 
(b) Inner part. 



This is compatible with an increase of Azo > 0.2 kpc in the 
case of i — 20° as described above. 



5.3 Heat increase 

We define the total heat in the disc as the kinetic energy of 
random motion 



^2 \ mi l Vl 



(14) 



where the sum goes over all particles in the disc, rrii is the 
mass of the ith particle, v< is the velocity of the ith parti- 
cle, and v c (Ri) is the mean circular velocity at the radial 
distance Ri of the ith particle. v c is first calculated in radial 
bins of equal particle numbers and then interpolated to the 
individual distances Ri. 
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Figure 10. (a) Total heat of the disc after merging (t = 1 Gyr) as 
a function of satellite's orbital inclination i for varying eccentric- 
ity e. -Eheat is measured in units of the heat of the isolated galaxy, 
Eheat ; so , at t = 1 Gyr. (b) Squared total velocity dispersion at 
t = 1 Gyr in the region around 8 kpc as a function of orbital incli- 
nation i and for varying eccentricity e. af ot is measured in units 
of the squared total velocity dispersion of the isolated galaxy. 



Figure [10] (a) shows Eh cat in units of the heat in the 
isolated disc £q lca t,iso at the end of the simulation (i.e. at t = 
1 Gyr). -Eheat, iso increases only by 0.18 per cent as compared 
to the initial value, demonstrating that our simulations are 
free of numerical heating. Figure \W\ (b) displays the same 
for the square of the total velocity dispersion in the solar 
neighbourhood. 

The heating efficiency of the satellite depends on its 
orbital inclination as well as on its orbital eccentricity. For 
low inclinations, satellites on less elliptical orbits heat more 
effectively than those on more elliptical ones. For i > 40°, 
however, satellites on orbits with e = 0.78 are more effective. 
Prograde satellites heat the most with an maximum increase 
of 22 per cent, while retrograde ones only have a small effect 
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Figure 11. Velocity dispersion ratio Oz/otj at t = 1 Gyr in the 
region around 8 kpc as function of orbital inclination i and for 
varying eccentricity e. 



of about 2 per cent. These results are consistent with those 
obtained bv I Velazquez fc White (19991 ). 

Observations show that the velocity dispersion in the so- 
lar neighbourhood increases from about 30 km s _1 to about 
60 k m s -1 in 8 Gyr approxim ately proportional to y/i (see, 
e.g. , lHanninen fc FIvnn 20021 ) . This means, that the specific 
heat increases linearly by about 1350 km 2 s~ 2 . If we assume 
that the heat increase is completely due to satellite mergers, 
then the required merger rate n is 



1 



1350 km 2 



" Ah " 8 s 2 Gyr ^ 

where Ah is the specific heat increase imparted by a single 
merger. 

In our simulations, we found that for low inclinations 
the heat increase is approximately uniform over the whole 
disc (except for the bulge-dominated inner part) and lies 
between 100 and 1000 km 2 s -2 per merger. For high incli- 
nations, the heat increase is mainly due to flaring of the 
outer parts of the disc (R > 15 kpc), while the solar neigh- 
bourhood is basically unaffected. Only considering low incli- 
nations, the required merger rate then lies between v = 0.17 
and v — 1.69 mergers per Gyr, depending on inclination and 
eccentricity of the orbit. 



5.4 Velocity dispersion ratio 

Observations show that the current ratio of the radial and 
vertical velocity dispersions in the solar neighbourhood, 
<y~//JR, is approxim ately 0.5 ±0.1 and incr eases slightly with 
stellar age, oc t ' 16 (|HoImberg et al. 2007ft . Figure [TTJ shows 
this ratio after merging as a function of orbital inclination i. 
The satellites with low prograde orbits decrease the ratio to 
at most about 0.4, i.e., they heat more efficiently in radial 
than in vertical direction. In one case [i — 20°, e = 0.56) 
there is an increase to approximately 0.6. For inclinations 
above 30° the ratio remains mostly unchanged in the solar 
neighbourhood. 



5.5 Energy and angular momentum transfer 

An interesting question is where the initial energy and an- 
gular momentum of the satellite end up in the final system. 

Figure [12] shows the change in energy of the various 
components. The energy is measured in units where the total 
energy of the isolated galaxy is —1. "Satellite" designates 
those particles that once made up the satellite, but now are 
distributed over the whole system. 

The initial galaxy-satellite system is not in equilibrium. 
As a consequence of virialisation, both the satellite and the 
halo gain kinetic and lose potential energy. This results in 
a deeper potential well in the centre, which reduces the po- 
tential energy of bulge and disc. During its in-spiral, the 
satellite transfers part of its kinetic energy to the halo par- 
ticles. This is due to dynamical friction. 

This effect can also be seen in Figure [TBI which shows 
the change in angular momentum as a function of inclina- 
tion, measured in units of the total angular momentum of 
the isolated galaxy. The angular momentum of the bulge 
barely changes. This is not surprising since the satellites are 
already mostly destroyed before reaching the inner part of 
the galaxy. Our data reveal that the radial distribution of 
angular momentum in the disc changes: the disc slows down 
and expands. Its total angular momentum, however, remains 
approximately nearly constant. It only increases slightly for 
satellites on prograde orbits and decreases slightly for satel- 
lites on retrograde orbits. The satellite's initial orbital an- 
gular momentum is 0.31 in these units. Accordingly, it loses 
between 15 per cent and 20 per cent of its angular momen- 
tum. Of that, about 80 per cent goes into the halo, while 
the rest is imparted onto the disc. 



5.6 Comparison with higher resolution 

To ensure that the choice of resolution does not influence 
the results we simulate a case with e = 0.56 and i = 0°, 
and higher resolution. We choose N = 256 and q = 0.125, 
as this is the highest possible resolution one can achieve on 
our present hardware. 

We calculate the Fourier transform of the density dis- 
tribution of the disc at the end of the simulation in order 
to find possible perturbations caused by low resolution. The 
mth complex Fourier coefficient is given by 



A m (Ri) = i-^Atf iC - im ^ 



where Ri and Si are the central radius and area of the 
ith bin, the index j runs over all particles in bin i, and Mj 
and ipj are the mass and polar angle of the jth particle. The 
amplitude of A m is denoted by A m and the phase by 9 m . 

Figure [14] shows the amplitudes of the first and second 
mode normalised by the mean density, while Figure [15] dis- 
plays the phase of the first mode. The two resolutions result 
in a roughly similar distribution in Fourier space. There are 
no strong perturbations in amplitude caused by the lower 
resolution. Neither is the pattern speed affected significantly. 



6 CONCLUSIONS 

Superbox is a particle-mesh code where additional grids 
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Figure 12. Change in total (a), kinetic (b) and potential (c) en- 
ergy of the four components with respect to the initial energies as 
a function of orbital inclination i. The energy change is measured 
in units of the total energy of the isolated galaxy. Only the data 
for e = 0.56 is shown. 
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Figure 13. Change in total angular momentum of the four com- 
ponents with respect to the initial angular momenta as a function 
of orbital inclination i. The change is measured in units of the to- 
tal angular momentum of the isolated galaxy. Only the data for 
e = 0.56 is shown. 
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Figure 14. Amplitudes of the first and second spiral mode as a 
function of radius at the end of the simulation. 



and sub-grids are applied to regions of high particle density. 
This strategy turns out to be very efficient since the code can 
run on any workstation or PC. Nevertheless, the code has its 
limitations. For instance, stellar discs are poorly resolved in 
vertical direction. We overcome this problem by introducing 
flattened grids. This is one of the features of the new code 
Superbox-10 where, in addition, an individual mass can be 
assigned to each particles. 

We found that the computationally most intensive part 
of the code is the FFT. We parallelised it using the library 
fftw. This resulted in a speed-up of 2 to 24, depending on 
grid size N and the number of particles n. 

We created a galaxy model with an exponential disc, a 
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Figure 15. Phase of the first mode as a function of radius at the 
end of the simulation. 



bulge with a Hernquist profile and a dark matter halo, also 
with a Hernquist profile. The model was realised using the 
program MaGaLie. 

The model was first simulated in isolation. These simu- 
lations show that flattening the intermediate grid is an effi- 
cient means to reduce numerical heating in the simulation. 

We also simulated the merging of the galaxy with small 
satellites in order to analyse the proposed disc heating due 
to the interaction. We find that satellites on prograde orbits 
with low eccentricity and inclination heat the disc most effi- 
ciently. If the heat increase in the solar neighbourhood were 
to be explained by that type of satellite mergers alone, a 
rate between 0.2 and 1.7 mergers per Gyr would be required. 
The detailed analysis of energy and angular momentum re- 
distribution shows that most of the satellites energy and 
angular momentum is transfered to the dark matter halo. 
This shows that the halo plays an important role even in 
50:1 mergers. This confirms, that simulations of such pro- 
cesses should represent the halo by live particles and not by 
a fixed background potential. The presented pilot study of 
high-resolution simulations of disc heating by merging satel- 
lite galaxies serves as a starting point for an extended pa- 
rameter study to quantify the heating rate in a cosmological 
context. 
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